function [pra, g, h] = profprob(a, DMES, EQ, V_R, lambda)

pra=zeros(length(a),1);
g=zeros(1,length(lambda));
h=zeros(length(lambda));

if size(DMES,1)>0
    
    spne = EQ(1,:)>=V_R;
    EQ(1,spne) = 0;
    EQ(:,~spne) = repmat([1;zeros(size(EQ,1)-1,1)],1,sum(~spne));

    peq = exp(DMES * lambda);
    peq = peq / sum(peq);
    
    % probability of observing each profile
    pra = EQ(a,:) * peq;
    
    if nargout>1

        % gradient of Pra
        g = EQ(a,:) * (peq .* (DMES - peq' * DMES));
    
        if nargout==3
            D = DMES - peq' * DMES;
            DD = (peq .* D)' * DMES;
            h =  - pra * DD + (EQ(a,:)' .* peq .* D)' * D;
        end
    
    end
    
end

